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The energy eigenstates and eigenvalues for the ID Grosss-Pitaevskii-equation for three different 
potentials have been determined numerically. For that a new and powerful algrorithm has been 
developed which is capable to work also in the regime of strong nonlinearities. The numerically 
found solutions completely agree with the case of known analytical solutions. 
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I. INTRODUCTION 



One of the most interesting problems in today's physics is the exploration of the quantum-gravity regime. This is 
due to the fact that General Relativity and quantum theory are not compatible which makes it necessary to search 
for a new theory called quantum gravity which at the end should lead to effective modifications of General Relativity 
and/or quantum theory. Another issue is that in some approaches gravity is regarded as a solution to the measurement 
problem in quantum theory. Therefore there are a lot of reasons showing that it is important to explore the interaction 
of quantum matter with gravity with better accuracy. One possibility to study the behavior of quantum matter in 
gravitational fields is neutron and atom interferometry One may even go further and investigate the energy 

eigenstates of quantum matter in a gravitational trap as has been pushed forward using ultracold neutrons at the 
ILL Q • In this experiment the various eigenstates manifest themselves through a neutron flux which depends on the 
height in a step-like form. One difficulty in this experiment is that the steps are of the order of /zm which comes from 
the strength of the gravitational acceleration. With the recently developed technology of Bose-Einstein condensates 
(BEC) in microgravity condition [1| another physical system is available for investigating the quantum-gravity regime 
for a wider range of parameters. Is is feasible to perform similar experiments with ultracold atoms in a Gravito- 
Optical Surface Trap (GOST) with a small and variable gravitational acceleration so that the density profile of the 
quantum states related to various energy levels can be measured with better resolution. 

The solution of the eingenvalue problem for the Schrodinger equation in such a GOST has been solved in terms of 
the Airy-functions in, e.g., [Q-Q- ^ n or der to be able to describe also the eigenstates for a BEC, we are solving here 
the eigenvalue problem for the nonlinear Schrodinger equations, that is, for the Gross-Pitaevskii equation (GPE). For 
doing so we developed in this paper a new numerical algorithm which is capable to find solutions of the GPE which 
belong to saddle points of the action. This algorithm is applied to three different potentials, the box, the harmonic 
. trap and the GOST. Our numerical solutions completely agree with the known analytic solutions for the box. For 
' a first description of the method and in order to present first results we restrict at the moment to one-dimensional 
! problems. 

In this paper in section [XT] we first state the problem and introduce the notation. In section IIIII we describe the 
newly developed algorithm and apply this method in section HVl for solving the energy eigenvalue problem of the GPE 
for three different physically relevant potentials. The paper closes in Section V with an outlook indicating further 
00 ' work in this direction. 
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II. THE MODEL 



We start with the time-dependent GPE which describes the dynamics of a BEC subject to two-particle interactions, 
given by the nonlinear term gsl^fe, i)| 2 , and to an external potential V CK t 

ihd t V( X ,t)= ^-^-A + V cxt (x,t)+g s \y( X ,t)\ 2 ^( X ,t), (1) 

where ^(x), x = (x,y,z) is normalized to the total number of particles TV = f n ^(x)! 2 ^ 3 ^. The GPE is valid for 
dilute condensates obeying the diluteness criterion, that is, the s-wave scattering length a and the average density of 
the gas h must fulfill h\a\ 3 <C 1. The nonlinearity parameter gs is determined by the scattering length via gs = 4 a , 
where m is the mass of the atom. Moreover, the scattering length can acquire both signs, having magnitudes of 
some nanometers. However, in this work we will focus on the case gs > which describes repulsive two-particle 
interactions. The function ^(x, t) has the meaning of an order parameter, is a classical field and is also interpreted 
as the wave function of the condensate. 
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For the calculation of the ground states and higher modes of a BEC in a time-independent external potential one 
makes the ansatz ^(x, t) = ^(x) exp (—ifj,t/h) leading to the stationary GPE 

p*(x)= f-|^A + ^ cx t(x)+.g5|*(x)| 2 ) *( X ), (2) 

where fi is the chemical potential. We also assume that the potential T4 x t(x) is bounded from below so that we can 
take F cxt (x) > 0. 

The stationary GPE can be derived from the action 

A[^;^]:=F[9]-^m, (3) 

with the free energy 



m ^ X ( V *W) 2 + ^ext(x)* 2 (x) + ^*<(x 



) d 3 X, (4) 



where we assumed a real \E'. The particle number is given by 

N[t>] := [ * 2 d 3 x. (5) 



Throughout this paper we will restrict ourselves to one-dimensional problems in order to demonstrate our new 
algorithm. 

In order to facilitate the numerical calculations, as one usually does, we rescale and renormalize the coordinates 
and the wave function according to 

x^Lx, *^%/iV*/i 3/2 , (6) 

where ^(x) is normalized to 1, leading to 

~ + V CK t(x)+ 1 ^ 2 (x)^j *(as) = e*(x), (7) 

with the dimensionless quantities 

- 2mL 2 V cxt (x) 2Nmg s 2m^L 2 

The length scale L is arbitrary and may depend on various physical parameters. It is chosen in such a way that the 
dimensionless quantities are convenient for numerical calculations. Note that in particular the nonlinearity parameter 
7 depends on the length scale L. 

Note that we do not restrict ourselves to functions that are normalized to one. Instead we are searching for solutions 
that are not normalized for a given pair 7, e. From equation ([7]) it is evident that each found solution can be normalized 
to one by adjusting the nonlinearity 7. 



III. THE ALGORITHM 
A. The general setting 

In computer numerics an attempt to solve nonlinear partial differential equations is to use some variant of the Newton 
method or the imaginary time propagation. The latter method is based on the the splitting and discretisation (e.g. 
Crank-Nicolson) of the unitary time evolution operator 

sua. This is reliable for ground state solutions. In this 

paper we will present a new Newton Method. 

Newton methods are gradient based algorithms that follow a descent direction until a local minimum of the action is 
reached. The direction depends on the choice of the inner product and/or of the preconditioning procedure. Solutions 
can be understood as critical points of some action A on the underlying dual space defined by the inner product. 
The type of a critical point related to a solution is determined by the eigenvalues of the Hessian [25[ evaluated at the 
critical point: 
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If all eigenvalues of the Hessian are positive then the critical point is a local minimum of A. 
On the other hand if all are negative then we have a local maximum. 

If we have a finite number of negative eigenvalues and all other eigenvalues are > then we have a horse saddle 
Fig. 1(b) In this situation the number of negative eigenvalues is the number of linear independent descent 



directions at this critical point. 

And the last case is when the Hessian is degenerated at a critical point, 
with a monkey saddle Fig. 1(a) for isolated critical points. 



For example this can be associated 



The number of negative eigenvalues is known as the Morse index. Solutions that belong to a local minimum of 
an action A are candidates for solutions which can be easily found by standard Newton methods, which searches in 
the whole L2 space. Unfortunately, finding critical points of a certain saddle type depends on an educated guess. In 
order to have a straightforward method at hand it is necessary to confine the search on a subspace of our Hilbcrt 
space. For linear eigenvalue problems this is easy to do because of the orthogonality of eigenfunctions. The gradient 
at every iteration step is orthogonalized with respect to the previously found eigenfunctions using the Gram-Schmidt 
procedure. Therefore it is easy to find eigenfunctions in ascending order of eigenvalues or, equivalently, in ascending 
order of the Morse index. Unfortunately, in the nonlinear case the orthogonality no longer holds, so that an other 
approach is needed. The basic idea is to constrain the quest for a solution to a submanifold in the underlying function 
space. 

In the following we need the first variational derivative (Gateaux derivative) 



A'[*;ii]h:= —A[V + eh;n} 



c=0 



which via 



(V La A[*;n],h) :=A'[*;p]h, 



(9) 



(10) 



can be identified with an L2 gradient V L 2 A[<j) k , /j]. Here (•, •) denotes the L2 scalar product. For the GPE we have 



(11) 



Therefore, if \& is a critical point of the action A then the L2 gradient of A vanishes and, hence, ^ is a solution of the 
GPE. 



B. Review of the Newton method 



The discrete Newton method is given by 

(j) k+1 = <t> k - Td k , 



(12) 
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where 

d k :=0- 1 V L2 A[(l) k ;fj,] (13) 

is the search direction and O^ 1 denotes a preconditioning operator that improves the convergence behaviour. In this 
context k is the iteration index and r the stepsize. The minus sign in front of r denotes that the correction of the 
step <p k is performed in the negative direction of the preconditioned L2 gradient. The stepsize can be a constant or 
can be determined at every iteration step using the linesearch or trusted region method. The solution then is given 
by cjf° l := lim^oe <f) k . In the nonlinear case the widely used Newton method is only capable to find Morse index 
zero solutions. Finding higher Morse index solutions for strong nonlinearities is a hard task and the standard Newton 
method is not able to do that. 

A von Neumann analysis applied to the standard Newton method (i.e. with the preconditioning O^ 1 = 1 ) leads 
to a convergence criterion like the famous Courant-Friedrich-Lewy condition which estimates a bound on the stepsize 
r depending on the discretization lengths and other parameters of the differential equation. Therefore a bad choice 
for r causes a failure of the Newton method. A too small r decreases the convergence rate. In order to handle this 
issue the preconditioning 0~ 1 is necessary. There are two well known methods, among others: 

1. A classical choice for O -1 is the inverse of the Hessian, or at least a numerical approximation. Due to the fact 
that computation time and storage space are precious and the full inverse Hessian is a dense matrix that is not 
fast computable new techniques have been invented to overcome this problem. The simplest one is to use the 
difference between two L2 gradients approximating the diagonal of the Hessian. 

2. A modern approach is to use the Sobolev preconditioning [Tl|, The L 2 gradient is mapped to a different 
Sobolev space, for example W 1 ' 2 . From a mathematical point of view the L2 gradient is filtered in Fourier space 
so that spatial oscillations are smoothed out. 

Upon the choice of preconditioning the direction of the gradient is altered. For a descend direction we have 
(0 _1 Vl 2 .A, Vl 2 „4) > for some arbitrary action A. The stepsize control of classical Newton methods fails if this 
condition does not hold. 



C. Newer approaches 

From equation (jj) it is evident that F[<fi k ] > for any <f) k 7^ and gs > so that <f> k = is the only critical point 
of F. Therefore it is only the term /j,N in ([3]) by which new critical points can appear. Accordingly, the key idea for 
the existence of solutions of non linear differential equations is to have terms that are capable of balancing the non 
linearity and all other terms. 

In order to emphasize the idea of balancing the nonlinearities we consider for demonstration purpose a classical 
case [T3I . [l4| for the situation without external potential and attractive two-particle interaction, thus V^xt — and 
7 < 0. Then the functional F reads 




F[0] = is a critical point and F is not bounded. Without further constraints a standard Newton method would 
fail. It is clear that there exists a t k 7^ that fulfils 

*W]tf fc = jf (^*Y jf(**) 4 <fc = 0, (15) 

which means that the kinetic part is balancing the interaction part. 

A Newton method which calculates the L2 gradient at the point t k (f> k defined by equation (Tl"5j) generates a sequence 
{t k ,(f> k } where the 4> k converge to a solution so1 7^ in L2, which is a local minimum of F. This is known as a 
minimization process restricted to the Nehari manifold 

t iei = extremum F[t k ^) k ] . (16) 

For the k-th step, t k ci cf> k is a reference point in the underlying function space where the L2 gradient is calculated. 
The definition (P) is equivalent to F' [t k cl (/) k }(f> k = 0. The sequence of functions (f> k calculated this way will converge 
to the solution so1 . With the restriction to the Nehari manifold it is possible to find Morse index one solutions of 
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([14]) . (For a general functional this may not always work.) Therefore it is, in general, very useful to confine the search 
for solutions to a manifold where the critical point lies in a local minimum on this manifold so that classical Newton 
methods are able to find such solutions. If one finds no extremum then the method is not applicable. 
A generalization of this idea has been presented in [l5j|. There, in the k-th iteration step the function 

D 

P? M ~^^T 4 +t fe fc , (17) 

i=l 

has been defined where t :— (t k . . .t k D ,t k j and the Tj are the previously found solutions, which were calculated by 
the algorithm [l5j], and D is the dimension of the support that is spanned by the T^. 

The idea behind this is to find solutions in the order of their Morse index which is similar to linear problems. First 
find the global minimum, then use the ground state to define a solution manifold in order to stay away from the 
ground state. After the the first excited state is found, it is used again together with the ground state to define a new 
solution manifold in order to stay away from the first and second solution. This is repeated until the algorithm fails. 
The ground state has Morse index zero and the first excited state has Morse index one. 

Then, for some action J[P?' k ] what can be interpreted as function of t, the extrema of the J[P?' k ] determine the 
vector t which is taken to define the reference point 

Cf = (t k ■ ■ ■ t k D ,t k ) = extremum J[pP' k ] . (18) 

Then the pP k is the reference point at which the L 2 gradient ^7l 2 <^[P^ > ' ] is calculated. 

If 7 > then the action may have minima. Zhou uses in [15j the so called active Lagrangian J := F— ^e k (N — 1) for 
his algorithm in order to find normalized solutions. The eigenvalue term e k is indispensable, but the final eigenvalue 
is not known from the beginning and has to be altered at every iteration step k. This includes the risk that at some 
iteration step k the solution is the trivial one t Ie f = 0. Zhou demonstrated that for 7 € 0(1) in a 2D GPE with a 2D 
harmonic trapping it is possible to find solutions in the order of their eigenvalues E{ where Si > £j_i. 

A related way to define a solution manifold [l6| is to require that the directional derivatives in direction of the D 
known solutions and the current iterations step vanishes. In order to find the reference point t TC f the following system 
of equations has to be solved: 

(v i2l 7[P^],T 1 ) = 
(v i2l 7[PP],T D ) = 

(v L2 ^[P^],0 fc ) =0. (19) 

This is a system of D + 1 equations for the D + 1 unknown variables (t k , t k D , t k ). The trivial solution is always a 
solution, but not the desired one. The numerical solution of this system depends on the initial guess of the vector 
(t k , t k D , t k ). Due to the nonlinearity of V L2^T[P?' k ) this system has more than one solution for suitable (e k ,<f) k ). 

In order to find an optimal reference point the initial vector (t k , tp, t k ) has to be guessed systematically. We used 
all corners and all center points of the faces of a D + 1 dimensional cube as guesses and took the initial guess where 
Vi<D\ti \ < \t°\. Under this condition we always succeeded to find a new solution for nonlinearities of the order 0(1). 
With increasing nonlinearity it becomes impossible to fulfil this condition when a solution with smaller non linearity 
is used as a guess. As a result of our numerical simulations, we like to add some remarks concerning the feasibility of 
the methods presented in [l5[ and [l6j : 

• First, for the reference point (|18|) and (fT9|) one encounters the problem that the old solutions also lie on the 
same manifold, so that if the initial guess is too far away from the final solution then one finds only old critical 
points. 

• Second, due to the fact that the eigenvalue e k is altered at every step there is a chance that it converges to an 
old one. 

• Third, in higher dimensions for non isotropic potentials the order of eigenvalues and eigenfunctions depending 
on the parameters of V ox t can be changed and the guess for the next solution may not be an appropriate one. 
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D. The new algorithm 

In order to overcome these problems, we developed a modified algorithm (see Fig. ((2J) which yields the following 
characteristics: 

1. Instead of using the previously found solutions "fi ■ • ■ "I^, we are calculating the reference point via 

(V L2 „4 [rfyn,. + M,.! A*n,.] . Vy s ) = 

(V L2 A [r 8 Vn,» + <£<,; M«, s ] , O - 0. (20) 

That means that we do not need the previous D solutions of lower Morse index as in Q19jl . Here r k , q k £ M. and 
fc and s are numerical counter variables which are only used in the algorithm. The quantum number n refers to 
the mode of the solutions we are interested in and is defined through the linear eigenfunctions. The functions 
ip n ,s and </>k s have the same nodal structure and 4>n s can be viewed as a correction to the solution ^„ iS for the 
previous eigenvalue. 

2. Unlike in the previously presented approaches we are working now with a fixed eigenvalue /i„ jS which is increased 
by a value A/i after a solution is found for the current fJ- n ,s, thus fi n .s+i = fJ-n.s + A/i. The nonlinearity 7„ iS is 
determined as a function of this chosen jU„ lS . 

3. The solution found in this way is not normalized to one. Therefore we normalize it and readjust the 7„ jS 
according to the particle number N. 

4. For the search direction d k , 

d k = 0~ l V L2 A [<.;/Xn,.] = (-^ + (Kxt - Mn, s ) + 3 7 „, s V i2 ^l M„, s ] , (21) 

we are using the inverse of the analytic Hessian evaluated at <^ s as the preconditioning operator O -1 . Our 
reference point defined by (f20|) together with the preconditioning pT]) assures that we do not leave the subspace 
with same nodal structure. In contrast, a Sobolev preconditioning would lead to ground state solutions only. 

The algorithm is implemented in CH — h The main part consists of two nested loops with the inner loop counter k 
and the outer loop counter s. The inner loop (STEP 1 to STEP 6) represents our Newton method. Within the outer 
loop /z„ iS is increased and storage operations are conducted. For convenience we take the calculated 7„ jS as starting 
point for the solution for the next eigenvalue fi„ tS +i. 

For the numerical derivatives a three point stencil is used. The integrals are evaluated with Simpson's rule and the 
differential equation from STEP 5 is solved in a finite difference setup with a Bi Conjugate Gradient solver. We used 
A/i = 0.5 , t = 0.01, dx = 0.05 and 1401 grid points. 

In the following we present the individual computational steps as depicted in Fig. @. 

STEP 1 The functions ipo an d 0° q are initialized to where ^ n is the analytic eigenfunction of the linear 
Schrodinger equation for the n-th quantum number. However, if ground state solutions are to be calculated, 
then one has to set ip n 0. In this case the solution manifold reduces to the Nehari manifold. At the end, set 

STEP 2 The numerical algorithm which solves the system of the two equations (|2U)) needs an initial guess for ££ ef . 
Due to the nonlinearity of Vl 2 .A s ; /x n ,J with respect to the s the solution is not unique. As guesses we 
used ^ cf = (0, 1) and tj? ef = (1,1) and selected for convenience the final ^ ef with the larger vector norm. 

STEP 3 Calculate the L 2 gradient using (ITTj) . 

STEP 4 Test for the quality of convergence. We use the maximum norm to check for convergence (the L 2 norm 
could be used also since both norms are equivalent). For our numerical calculations we set r\ = 1 • 10~ J . On 
expense of more iteration steps better results can be achieved for smaller r\. 

STEP 5 In general, the numerical inversion of the operator O on the l.h.s. of (l2Tj) consumes much computer storage 
and calculation time. In order to avoid this, we solve the differential equation of STEP 5 in figure [5] and, thus, 
obtain the search direction d. The disadvantage of this procedure is that the differential operator O has to be 
assembled at every iteration step. 
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STEP 1 

Calculate the Eigcnfunction for 7 — 0. 
Calculate the Eigenvalue E n for 7 — 0. 
Set s <- 0, 0° <- *„, Vn,o <- *n 



STEP 2 

if fe — then find optimal ref. point i^ ef else 
find the new ref point using i^T 1 as a guess. 



STEP 3 



Calculate V L2 A[Pl ,k ] 



STEP 4 



STEP 5 

Solve: 



fa2 + (V cxt - Mn.a) + 3 7 „,a(P3'* ) 2 1 d fc 




STEP 6 

= 0^_ s - T sgn(q k )d k 
k <- k + 1 



yes 



STEP 7 

Calculate ./V 
Store (/i„, 3 ,7„, s iV, jV _1/2 Pi' fc ' 
7n,s + l <- 7n,»JV 
Replace ?/j n s with _Pl ,fc 

*rcf 

s «- a + 1 , fc <- 

if < M„,f in 



yes 
STEP 8 

Exit 



FIG. 2: Flow chart of the algorithm 



STEP 6 A modified Newton step is carried out. We set the stepsize to a constant value r = 0.01. A classical linesearch 
or trusted region stepsize control is not applicable due the fact that d is not always a descent direction. Note 
that equation ([20)1 is invariant under the simultaneous sign change of r k and q k . This reflects the invariance of 
the GPE under the transformation of the wave function i/j — > —ip. Thus, the factor sgn(q^ ) has to be introduced 
into the second term in the r.h.s. of (1121) in order to have a unique notion of ascent and descent directions, 
respectively. Therefore, the search direction can be made unique by means of multiplying r<i[<^ s ; fi nyS ] with 
sgn(gj). 

STEP 7 First we calculate the particle number N of Pi' k according to ([5]). Then we save the solution Pl' k for the 

current eigenvalue n n ,s, the adjusted nonlinearity 7 S , and the corresponding normalized solution N~ 1 / 2 Pi' k . 

After that we replace the function -i/>n,s by the just calculated P~' k - Before we proceed to the next inner loop 
s + 1 we increase the eigenvalue ^„ jS by A/x, where A/i has a typical value of 0.5. Go to STEP 2. 



In Figs. 3(a)|3(c) we show typical forms of the error estimate ||Vl 2 -4[<^ s ; /i„ lS ]||oo as a function of the inner loop 



counter k for each of the three problems discussed later. The number of iteration steps for this algorithm applied 
to these problems was of the order of O(10 2 ) — O(10 3 ) for the inner loop. From these graphs it is evident that the 
error estimate is not necessarily decreasing right from the first iteration step k = as one might expect. Thus our 
algorithm permits that the search direction (Vl 2 -4 s ; n n ,s] > d k ) can be descending or ascending in contrast to 
the aforementioned algorithms (see Appendix IA1I) . We have to emphasize that this depends on the preconditioning. 
In these logarithmic plots the linear behaviour reveals the exponential decay of the norm. 



IV. SOLUTIONS FOR VARIOUS POTENTIALS 



In this section we present analytical and numerical solutions for the energy eigenstates and the energy eigenvalues 
for the GPE for three potentials, that is, (i) for a box, (ii) gravitational surface trap, and (iii) the harmonic trap. While 
usually in experiments BECs are created in the ground state, excited states might emerge through an appropriate 
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(a)n = 5 for the box. 




10 15 20 25 30 35 40 
k 

(b)n = 5 for the COST. 



3. 



> 




1 2 3 4 5 6 
k 

(c)n = 5 for the harmonic trap. 

FIG. 3: Maximum norm of the L2 gradient as a function of the iteration counter k for the first five outer loops iterations. 



periodic motion of, e.g., the walls of a box potential. This is similar to the creation of waves of a viscous fluid in a 
box through the motion of walls. The explicit procedure of the creation of excited states of a BEC obeying the GPE 
will be discussed in a subsequent paper. 



A. BEC in a box 

In this section we present the numerical results of a BEC confined in a box of finite size L. With ([6]) and the natural 
length scale L = h/ y/2mfl the dimensionless GPE reads 

("* 

The potential is given by 



-^2 + tUfc) + 7*nW ) *n(x) = en*n(x) . (22) 



Vext{X) = < . (23) 

I 00 else . 

As usually, we require the standard boundary conditions \l/ n (0) = and *n(l) = 0. 
For 7 = the eigenfunctions and energies are simply given by 

^n( x ) = V~2 sin(imx) and e n = n 2 n 2 , (24) 

where n = 1, 2, 3, . . .. 
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For 7 > this problem can be solved analytically by means of the Jacobi elliptic function sn 



^n(x) = 2\/2m7 1 nK(m)sn (2nK(m)x\m) 



(25) 



where K is the complete elliptic integral of the first kind. The definitions of the elliptic integrals and functions are 
taken from [l8j]. 2nK(m) is the real period of the Jacobi sn function. The modulus m of the Jacobi sn function is 
determined by the following equation 



8n 2 (K(m) - E(m)) =7, 



(26) 



which is derived from the normalization condition and E(m) is the complete elliptic integral of the second kind. The 
energy spectrum then is 



4n 2 K(m) 2 (l + m). 



(27) 



The limiting case 7 = leads to m — 0, K(0) — ^, E(0) = ^, and the Jacobi elliptic function sn in equation (f2"5)l 
reduces to sin(7rna;). 




3, 




0.2 0.4 0.6 0.8 



(a)/x = 23, 7 « 9.1865 



(b)/i = 54, 7 « 9.8271 



FIG. 4: Comparison of numerical solution (solid red line) and the analytic solution (dots). 



The knowledge of these analytically given solutions is very useful as a benchmark for our algorithm. In Figs. 4(a) 
and |4(b)| the comparison between the analytical and numerical solution for the ground state and the first mode is 
shown. The solid red line is the numerical result and the dots are calculated with the analytic solution. The deviations 
are of the order of O(10~ 4 ). For large nonlinearities the numerical calculation of modulus with equation (f2"6")l becomes 
difficult because the number of required decimals increases fast. 

In Fig. [5] the first six modes are plotted, starting with the strictly positive zeroth mode. The solid black lines show 
the cigcnfunctions for the linear case 7 = 0. The other two lines show the numerical solutions of the nonlinear problem 
for different eigenvalues All solutions are normalized to one. The corresponding eigenvalues fi n for a given 7„ can 
be read off from table U or from Fig. [6] The relation between 7„ and /i„ in Fig. [6] is proportional but not linear. 

The amplitudes of the wave function shown in Fig. [5] decrease with increasing j n and the maxima and minima 
become more and more flat. This can be easily understood from the fact that the repulsion becomes stronger for larger 
7„ so that the wave functions tend to a spatial equalization. As a consequence, the gradient of the wave functions at 
the boundary grows and with that the kinetic energy. 

Note that for increasing mode numbers and at fixed nonlinearity 7 the broadening effect gets smaller. This can be 
seen from the solutions for /io = 500 at mode zero and for the 5-th mode at ^5 = 1000 (see Fig. [5]). 



B. Gravitational Trap 

Now we solve the GPE with a gravitational potential. With the potential 

T _ , . I max if x > 
Kxt (x) = { 

00 else . 



(28) 



10 



mode 



mode 1 



1 


6 


1— ' 


4 


(— * 


2 




1 





8 





6 





4 





2 










2 







1000 






500 











0.2 0.4 0.6 0.; 

x 





mode 3 




mode 4 



mode 5 





FIG. 5: The first six solutions of the GPE in a box. The corresponding nonlinearities j n for given fj, n can be found in table (JT|. 



we have the natural length scale L = (h 2 /2m 2 g) so that the dimensionless GPE reads 

(~^2 + x + ^n(x) 2 ^J V n {x) = e n ^ n (x) . (29) 

In the linear case with 7 = the eigenfunctions are well known Q ■ The general solution is a linear combination of 
the AiryAi and AiryBi functions where the AiryBi is omitted since it is not compatible with the boundary conditions 
x I / n(0) = and ^ n (oo) = 0. The n-th eigenfunction is given by 

*„0) = A n Ai(x + x n ) , (30) 

where x n is the n-th zero of the AiryAi function and of course the orthogonality relation (ty n (x), ^ m (x)} = S nm holds. 
All zeros are negative so that the normalizable part of the general solution is shifted to the right. The n-th eigenvalue 
e„ is also given by the n-th zero. Unfortunately no analytic expression for the normalization factor A n exists. Hence, 
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FIG. 6: The one-to-one correspondence between the energy eigenvalues fi and the nonlinearity parameter 7 for the box potential 
for different eigenmodes. The solid line represents the ground state, increasing mode numbers to the right. 



TABLE I: 7„(ai„) for BOX 



mode Jnjpm — 500) 7n(^n = 1000) 

436.7686 910.5884 

1 373.5441 821.2079 

2 310.1694 731.8231 

3 245.2187 642.2824 

4 175.4825 551.4566 

5 98.0917 457.6764 



it is given by 



dxAi(x + x n ) 2 . (31) 

In the nonlinear case it is much more complicated to find solutions because equation (|29p admits a huge number of 
not normalizable solutions. Most of them have poles on the real axis. 

For e n = and 7 = 2 problem (|29|) is known as the Painleve II equation (PII). This equation possesses the Painleve 
property [19j which is a condition of integrability. These differential equations cannot be integrated by means of 
elementary functions. A possibility of finding solutions to PII is to solve the corresponding Riemann-Hilbert problem 
numerically [2(j ■ Using this method many different solutions can be found, including nonphysical ones [Hj]. 

In Fig. [7] the first eight solutions of the GPE with linear potential are calculated for different values of /i„ with our 
new method. The solid black lines depict the AiryAi solutions for the linear case. The corresponding nonlinearity 
factors can be found in table [TTJ The difference here is that the mathematical domain is not finite and that for x — > 00 
the potential is diverging. However, in the numerical implementation the domain has finite size, i.e. of length L. On 
contrary to the standard treatment of boundary conditions it is only necessary to specify the value at ^(0) due to 




TABLE II: 7„(/x n ) for V4xt = x 



mode 


7n(Mn = 10) 


7«(m« = 20) 


7n(/ttn = 30) 





45.4993 


193.6644 


442.2494 


1 


36.8594 


181.1653 


426.8558 


2 


28.5681 


168.8259 


411.5684 


3 


20.6737 


156.6578 


396.3996 


4 


13.2327 


144.6551 


381.3278 


5 


6.2966 


132.8302 


366.3652 


6 




121.1916 


351.5141 


7 




109.8536 


336.4804 
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FIG. 7: The first eight solutions of the GPE for the GOST. The corresponding nonlinearities 7„ for given /i n can be found in 
table (Hill. 
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the diverging nature of the potential the value at *&(L) adjusts itself automatically. We do not pose any boundary 
conditions explicitly for ^(L) so that the finite size domain has no effect on the solutions. 

The first observation is that for a given mode number n > higher nonlinearities cause a quenching of the region 
between the boundary and the outermost maximum in comparison with the linear case. This can be understood by 
taking into account that V(x = 0) = oo limits the space on the left side for encountering particles. As a result, they 
move towards the outermost maximum causing a depletion of the particle number on the left. 
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FIG. 8: The one-to-one correspondence between the energy eigenvalues /j, a and the nonlinearity parameter 7 for the GOST 
potential for different eigenmodes. The solid line represents the ground state, increasing mode numbers to the right. 

The second observation which seems to be surprising is that the bulk of the wave function for different modes 
appears not to be changing for fixed fi n . However the explanation is simple: with increasing modes the nonlinearity 
parameter 7„ decreases at fixed [i n . This behaviour can be clearly seen in Fig. [5] Higher nonlinearities always enlarge 
the bulk of the solution. Fig 7„(/i„) in Fig. [8] shows the one to one correspondence between 7 and fi. 



C. Harmonic Trap 



As a last example we discuss the BEC in a harmonic potential given by 

1 



V cxt = ^u 2 x 2 . (32) 



With the natural length scale L = i/fi/mu the dimensionless equation reads 

d 2 



it 2 +x 2 + j^ n (x) 2 ) tf n (a;) = e„*„(a;). (33) 



In the linear case 7 = the solutions are 



*„(a;) = . = ^ exp (-.t 2 /2) H n (x) , (34) 



where 

(ki- 



rn 

H n (x) = (-1)" exp (x 2 ) jL. exp (-x 2 ) (35) 



are the weighted Hermite polynomes. 

As far as we know there are no analytic solutions known for this potential for 7 7^ 0. Numerical simulations basically 
focus on the zeroth mode 0, [22j ■ For the zeroth mode there is a rough approximation which can be obtained by 
neglecting the kinetic energy term in equation (|33j) . Then we have an algebraic equation which can easily be solved 
for ^oi x ) an d is known as the Thomas-Fermi solution 
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In Fig. [TU] we compared the Thomas-Fermi solution with the numerical solution of the ground state. For large 
nonlinearities the Thomas-Fermi approximation agrees very well with the numerical results in the center region of 
the condensate. 
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FIG. 9: The one-to-one correspondence between the energy eigenvalues fi„ and the nonlinearity parameter 7„ for the trapping 
potential for different eigenmodes. The solid line represents the ground state, increasing mode numbers to the right. 
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FIG. 10: Comparison between the num. solution (solid red line) and the Thomas Fermi approximation (blue dashed line). 

In Fig.[TT]thc first eight numerical solutions for the harmonic oscillator potential are depicted for different eigenvalues 
jj, n calculated with our new method. The corresponding nonlinearities 7„ can be found in table (IIIII) . The solid black 
lines correspond to the linear case with the weighted Hermite polynomial. In the numerical implementation there are 
no boundary conditions specified. The diverging nature of the potential forces the wave function to decay for x — > 00. 

The curves in Fig. [TTJ show that with increasing nonlinearity particles from the center region are pushed towards 
the outer region whereas the inner structures are squeezed. The harmonic potential has a much higher confinement 
so that the bulk remains relativity small compared to the gravitational potential. 



V. DISCUSSION AND OUTLOOK 



In this article we presented a new algorithm that is capable to find higher Morse index solutions of the stationary 
GPE for large nonlinearities. Mathematically speaking, these are saddle point solutions. The three crucial points 
are (i) to start with a fixed eigenvalue fj,, (ii) the reference point that contains only a function of the same Morse 
index in the support and (iii) the preconditioning of the Li gradient by using the analytic expression for the Hessian. 
Furthermore we demonstrated that we can find solutions for the GPE with large nonlinearity parameter in external 
potentials by starting only with the eigenfunction for the n-th mode of the corresponding linear problem. 
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FIG. 11: The first eight solutions of the GPE with the harmonic trap potential. The corresponding nonlinearities 7„ for given 
\i n can be found in table 
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TABLE III: j n (fj, n } for T4 xt = x 2 



mode 7nQ<n = 10) 7n (/fn = 20) 7n (M'l = 30) 7n(Mn = 40) 7n(Mn = 50) 

41.6008 118.0873 218.6985 336.9562 471.0779 

1 32.5650 106.1331 203.1434 319.0604 451.0741 

2 23.4540 93.4226 187.6607 301.1403 431.0489 

3 14.2610 80.7247 172.1164 283.2037 411.0061 

4 4.8558 68.0688 156.5821 265.2629 390.9532 

5 55.5007 141.0820 247.3326 370.8997 

6 43.0556 125.6424 229.4285 350.8563 

7 30.7438 110.2917 211.5678 330.8343 



In summary we calculated the eigenfunctions and energies for the one dimensional GPE for three classical potentials: 
the box, the harmonic trap and the GOST. In the case of the GOST we obtained higher order modes up to order 
seven for large nonlinearity parameters in the range of 7 = 336 — 442. To our present knowledge this seems to be 
the first time that solutions to the GOST setup for such high nonlinearities and high modes have been calculated. 
Furthermore, we showed that in the case of the box the numerically found solutions completely agree with known 
analytical solutions which confirms our algorithm. Also, there is a good agreement between the numerical zero mode 
solutions for the trap and the Thomas- Fermi approximation in the central region of the BEC. 

The next logical step is to apply this algorithm in 2D and 3D setups of the aforementioned three cases so that 
more realistic physical systems will be modelled. Moreover, the physical stability may be checked by propagating 
the solutions in time. Another issue which may be treated in future is to include self-gravity effects. At first, for 
an efective equation as the GPE self-gravity should be considered. For very dilute gases one may expect no effects 
but for high density BECs corresponding effects should be estimated. Self gravity also is an idea stated by Penrose 
[23| to understand the collapse of the wave function. Therefore it might be of interest to investigate whether in this 
context such effects might be accessible to experiment. We also plan to adopt our method to the case of coupled 
many component GPEs. 

Finally, concerning the algorithm, in the future it may be interesting to find a new stepsize control that incorporates 
the ascent direction. This may improve convergence behaviour. For spatial dimensions larger than 1 we also would 
like to extend our algorithm in order to incorporate also different coordinate systems which are more adapted to the 
physical problem. 
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Appendix A: Search direction 

For the action A we can write down the Taylor expansion around (f> k up to first order: 

A ; Ms ] = A [0*; Ms ] + A' {4> k+1 -4> k )+0 ({4> k+1 - <t> k f) . (Al) 

Using the Newton step (fT2"]l and the search direction (f2"Tj) equation (|A1|) can be written as: 

A[cp k+1 -^ s ] =A[4> h >,n a ]-TA! [0 k - lf i s ] d k + 0(d 2 ) 

= A[4> k -^s\ ~r{V L2 A[4> k -^s] ,d k )+0{d 2 ). (A2) 

If A [c/) k+1 ; fi s ] < A[<j) k ;fi s ] and (V l 2 A [<j) k ; fi s ] ,d k ) > then ~d k is a descent direction. If A [(j) k+1 \ fi s ] > 
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A[cj) k -ii s ] and (V L2 A[(/) k ;n s ] ,d k ) < then -d k is a ascent direction. 
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